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A recent paper by Lucas-Serrano et al. [1] indicates that a high-resolution central (HRC) scheme 
is robust enough to yield accurate hydrodynamical simulations of special relativistic flows in the 
presence of ultrarelativistic speeds and strong shock waves. In this paper we apply this scheme 
in full general relativity (involving dynamical spacetimes), and assess its suitability by performing 
test simulations for oscillations of rapidly rotating neutron stars and merger of binary neutron 
stars. It is demonstrated that this HRC scheme can yield results as accurate as those by the so- 
called high-resolution shock-capturing (HRSC) schemes based upon Riemann solvers. Furthermore, 
the adopted HRC scheme has increased computational efficiency as it avoids the costly solution 
of Riemann problems and has practical advantages in the modeling of neutron star spacetimes. 
Namely, it allows simulations with stiff equations of state by successfully dealing with very low- 
density unphysical atmospheres. These facts not only suggest that such a HRC scheme may be a 
desirable tool for hydrodynamical simulations in general relativity, but also open the possibility to 
perform accurate magnetohydrodynamical simulations in curved dynamic spacetimes. 
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Hydrodynamics simulation in general relativity (GR) 
is the best theoretical approach for investigating dynam- 
ical phenomena in relativistic astrophysics such as stellar 
core collapse to a neutron star and a black hole, and the 
merger of binary neutron stars. In the past several years 
this field has witnessed major development, to the stage 
that it is now feasible to perform accurate simulations 
for such general relativistic phenomena (see e.g. [2-6]). 
Currently, the most favored approach to hydrodynamics 
simulations in full GR combines the use of the so-called 
BSSN formalism to solve Einstein's field equations [7] and 
upwind high-resolution shock-capturing (HRSC) schemes 
to solve the hydrodynamics equations [8] in conserva- 
tion form. Hereafter, HRSC schemes are referred to as 
those in which the hydrodynamics equations are solved 
by means of (either exact or approximate) Riemann 
solvers [9,8] (i.e. Godunov-type schemes). 

Regarding the solution of the hydrodynamics equations 
it has been shown in a few recent papers [10,1] that high- 
resolution central symmetric schemes (HRC scheme here- 
after) yield numerical solutions as accurate as those by 
HRSC schemes for special relativistic flows (see e.g. [11] 
for a general introduction to HRSC and HRC schemes). 
The main conclusion of those works highlights the im- 
portance of the conservation form of the adopted scheme 
(either upwind or central) in conjunction with high-order 
cell-reconstruction procedures (to compute the numeri- 
cal hydrodynamical fluxes at cell interfaces) to gain ac- 
curacy while reducing as much as possible the inherent 
diffusion of central schemes at discontinuities. It is well- 
known that if a numerical scheme written in conserva- 



tion form converges, it automatically guarantees the cor- 
rect Rankine-Hugoniot (jump) conditions across discon- 
tinuities. This shock-capturing property is hence shared 
by both upwind and symmetric schemes. For practical 
reasons the most appealing feature of HRC schemes is 
the fact that, contrary to upwind HRSC schemes, they 
entirely sidestep the use of Riemann solvers, which re- 
sults in a great simplification for their numerical im- 
plementation as well as in enhanced computational ef- 
ficiency. However, it has not yet been clarified whether 
HRC schemes can also yield numerical results as accu- 
rate as those of HRSC schemes for simulations in full 
GR involving dynamical spacetimes. 

The aim of this paper is to demonstrate the robust- 
ness of a particular HRC scheme proposed by [12], and 
first used in special relativistic hydrodynamics by [1], for 
problems in full GR. As we have done in previous pa- 
pers (e.g. [2,3,13]), test simulations in both axisymmetry 
(rotating neutron stars) and full three-dimension (binary 
neutron star mergers) are performed to assess this fact. 

The numerical simulations are carried out using the 
same mathematical formulation as in [5], to which the 
interested reader is addressed for details about the basic 
equations, the gauge conditions, and the computational 
method. Einstein's evolution equations are solved using 
the so-called BSSN formalism [7] , adopting a slight vari- 
ation of the original form of the equations, which is re- 
ported in [5]. The hydrodynamics equations are written 
in conservation form and solved using both a Roe-type 
HRSC scheme [13] and a HRC scheme [1], with either 
the PPM third-order cell-reconstruction or the MC slope 
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limiter. Violations of the Hamiltonian constraint and 
conservation of ADM mass and angular momentum are 
monitored to check the accuracy of the simulations. 

We use a fixed uniform grid for both the axisymmetric 
and the three-dimensional (3D) simulations. The former 
are carried out in cylindrical coordinates (137, z) assum- 
ing equatorial plane symmetry. Computational grids of 
size {N + 1,N + 1) with N = 90, 120, 180, 240, and 
360 are used, with which convergence is shown. The 3D 
simulations are performed in Cartesian coordinates as- 
suming equatorial plane symmetry as well. In this case 
the grid adopted in the present test simulations consists 
of (377,377,189) zones for {x,y,z) respectively. 

In the axisymmetric simulations of isolated rotating 
neutron stars a F-law equation of state (EOS) is used, 
i.e. P = {T — l)pe. Here, P is the pressure, p the rest- 
mass density, e the specific internal energy, and F the 
adiabatic constant for which we choose the values 2 and 
2.5. The initial conditions for the equilibrium models are 
built using a polytropic EOS P = Kp^ , where K is the 
polytropic constant. 

Correspondingly, for the 3D simulations of binary neu- 
tron star merger a hybrid EOS is adopted, as described 
in [6]. In this EOS, the pressure and the specific inter- 
nal energy are written in the form P — Pcoid(/o) + Pth 
and £ = ecoid(/3) + £th where Pcoid and Ecoid are the cold 
(zero-temperature) parts, and are functions of p only. 
On the other hand Pth and £th are the thermal (finite- 
temperature) parts. During the simulation, p and e are 
evolved, and thus eth is determined by £ — £coid- For Pth> 
we simply set Pth = (Eth - l)p£th with Pth = 2. For 
the cold part of the hybrid EOS we use realistic EOS for 
zero-temperature nuclear matter, more precisely the SLy 
EOS [14]. 

As customary in grid-based hydrodynamics codes an 
artificial low-density atmosphere needs to be used in 
those regions outside the star representing vacuum. The 
density has to be low enough so that its presence does 
not affect the actual dynamics of the star. In previous 
simulations using a Roe-type HRSC scheme, a uniform 
density atmosphere as low as patm = 10~^/3max was used, 
where Pmax is the maximum density. (For soft EOS this 
value can be much smaller; e.g., Patm < 10~^^/9max for 
F = 4/3.) Lower values for the density could result in 
numerical instabilities developing around the stellar sur- 
face. However, we have found that when using the HRC 
scheme the threshold density in the atmosphere can be 
much smaller. The results presented next for the HRC 
scheme correspond to patm = 10~^°/9max, irrespective of 
the EOS used. 

We start discussing axisymmetric simulations of oscil- 
lations of rotating neutron stars. For these simulations 
we build rapidly rotating neutron stars with uniform an- 
gular velocity. This velocity is chosen so that it reaches 
the Kepler (mass-shedding) limit at the equatorial stellar 
surface. 

Two rotating neutron star models are considered. In 
one case F = 2 and the baryon rest mass is 90% of the 




FIG. 1. Evolution of the central density in units of the ini- 
tial value for the two rotating neutron star models considered. 
Time is shown in units of the rotational period of the neutron 
stars P. The dotted-dashed, dashed, dotted, and solid curves 
denote the results with (121,121), (181,181), (241,241), and 
(361,361) grid resolutions, respectively. 



maximum allowed value for uniformly rotating neutron 
stars of identical EOS. This model is the same as model 
R2 in Ref. [13], which allows for a direct comparison. The 
other model corresponds to F = 2.5 and is 95% of the 
maximum allowed mass. This is a very compact model, 
since the compactness parameter, defined as GM/Rc^ 
where M and R are the ADM mass and circumferential 
radius around the equatorial surface, is 0.214. For both 
models, the axis ratio of polar radius to equatorial radius 
is about 0.6. The ratio of the coordinate radius of the 
outer boundary of the computational grid to the stellar 
coordinate radius at the equator is 3. The simulations 
are started by reducing the pressure by 1% uniformly. 

Figure 1 shows the time-evolution of the central density 
for these two models obtained using the HRC scheme for 
the hydrodynamics equations. Each curve corresponds 
to a different grid resolution as explained in the caption. 
It is found that the HRC scheme succeeds in keeping the 
stars in equilibrium in such a dynamical spacetime. The 
neutron star oscillations can be followed accurately for 
more than 20 rotation periods. With small grid sizes 
(dotted and dashed lines) , the density experiences a sec- 
ular drift, decreasing with time gradually. The reason is 
that the angular momentum of the star is transported 
outward by numerical diffusion. However, this drift de- 
creases with improved grid resolution, and with the high- 
est resolution the average value of the central density is 
kept approximately constant. Second-order convergence 
is also achieved. 

It is worth to emphasize that despite the use of an arti- 
ficial atmosphere of tiny density, the HRC scheme makes 
it possible to follow the evolution of compact neutron 
stars with stiff EOS with F = 2.5 and to compute their 
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FIG. 2. Evolution of the ADM mass (top), angular momen- 
tum (middle), and violation of the Hamiltonian constraint 
(bottom) for F = 2.5. Time is shown in units of the rota- 
tional period of the neutron stars. The dashed, dotted, and 
solid curves denote the results with (91,91), (121,121), and 
(181,181) grid resolutions, respectively. 



fundamental oscillation frequency. Such simulation has 
not yet been accurately performed with HRSC schemes. 

In Fig. 2, we show the evolution of the ADM mass, 
angular momentum, and the averaged violation of the 
Hamiltonian constraint (in which the baryon rest mass 
density is used for the weight; see [5] for definition) for 
r = 2.5. (Similar results are obtained for F = 2.) The 
figure shows that the conserved quantities remain con- 
served to high accuracy, particularly for the finest grid, 
and that the violation of the Hamiltonian constraint re- 
mains small. The outstanding feature is that the depar- 
ture from angular momentum conservation with the HRC 
scheme is much smaller than with the HRSC scheme (see 
e.g. Fig. 5 in [13]). In the previous implementation the 
angular momentum gradually increases with time mainly 
due to the numerical error generated around the stellar 
surface and in the low-density atmosphere for which an 
artificial friction term was added to stabilize the compu- 
tation. With the HRC scheme, such a drift in the an- 
gular momentum conservation is suppressed within 0.1% 
error after 20 rotation periods for grid resolutions with 
N > 180, probably due to less numerical inaccuracies 
around the stellar surface. These results indicate that 
the HRC scheme used is a robust scheme for the simula- 
tion of isolated neutron stars. 

We now turn to present the results of numerical sim- 
ulations of binary neutron star mergers. In the present 
test we choose binaries of equal mass with 1.3-1.3Mo 
and 1.4-1.4Mo. As found in [6] using a HRSC scheme, a 
massive neutron star and a black hole are formed for the 
former and latter cases, respectively. 

In Fig. 3, we show the time-evolution of the central 
value of the lapse function, ac, and the maximum value 
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FIG. 3. Evolution of the maximum value of p and the cen- 
tral value of a for the merger of equal-mass binary neutron 
stars, 1.3-1.3M0 and 1.4-1.4M0. In the former and lat- 
ter cases, a massive neutron star and a black hole are the 
end-products of the merger, respectively. 



of the density, Pmax, for these two models. The solid and 
dashed curves indicate the results obtained with the HRC 
and HRSC schemes, respectively. In the smaller mass 
case, a massive neutron star is formed after the merger, 
and hence, ac and Pmax show a series of small-amplitude 
oscillations until they eventually relax to quasi-stationary 
values. For both hydrodynamical schemes the amplitude 
and the frequency of the resulting neutron star oscilla- 
tions agree well with each other. On the other hand, 
the outcome of the merger of the 1.4-1.4M0 binary is 
a black hole, as can be directly inferred from the rapid 
collapse of the central lapse and the rapid growth of the 
maximum density (from t ^ 2 ms onwards). Black hole 
formation is signaled by the appearance of an apparent 
horizon, which is detected in both implementations. In 
particular the time of formation of the apparent horizon 
agrees approximately for both schemes, with a time dif- 
ference of about 0.07 ms. For the two binary mergers 
considered, a small time lag in the evolution of ac and 
Pmax is observed between the two results computed by 
the different schemes. Its origin is likely the difference 
in the magnitude of the friction term around the stellar 
surface already discussed before which could generate an 
error in the angular momentum conservation. As men- 
tioned above, this error is smaller with the HRC scheme. 

In Fig. 4, we show the evolution of the ADM mass, an- 
gular momentum, and averaged violation of the Hamil- 
tonian constraint for the 1.3-1.3M0 binary merger. In 
this case, the ADM mass and angular momentum of the 
system show a gradual decrease due to the emission of 
gravitational waves [6]. Again, it is found that the two 
results agree well with each other within ~ 0.5%. The 
averaged violation of the Hamiltonian constraint remains 
approximately of identical magnitude, ~ 0.02, which in- 
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ent when the wave structure of the hyperbohc system to 
solve is unknown, as it is partiaUy the case in general rela- 
tivistic magnetohydrodynamics (GRMHD). Hence, HRC 
schemes can help the achievement of GRMHD simula- 
tions in which the equations to solve are more compli- 
cated than those of purely hydrodynamical flows [15]. 
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FIG. 4. Evolution of the ADM mass (top), angular mo- 
mentum (middle) , and violation of the Hamiltonian constraint 
(bottom) for the merger of binary neutron stars with masses 
1.3-1.3Mq. The solid and dashed curves in all panels indi- 
cate the results obtained with the HRC scheme and with the 
HRSC scheme, respectively. 

dicates that the accuracy of the results of the two hydro- 
dynamical schemes is approximately identical. 

To summarize, it has been shown through simulations 
of pulsating and rotating neutron stars, and binary neu- 
tron star mergers, that the results produced by the HRC 
scheme proposed by [12] agree well with those obtained 
with a Roe-type HRSC scheme. The accuracy measured 
by the evolution of the ADM mass, angular momentum, 
and violation of the Hamiltonian constraint in the HRC 
scheme are as good as or even better than those obtained 
for the HRSC scheme. In addition, the HRC scheme 
has a number of advantages to the HRSC scheme: (1) 
it is straightforward to implement since the solution of 
Riemann problems is avoided; hence one does not need 
to compute the complicated sets of eigenvectors of the 
Jacobian matrices associated with the fluxes (transport 
terms) of the hydrodynamics equations; (2) for this rea- 
son the computational costs of the HRC scheme are much 
less expensive, as the characteristic information required 
in HRSC schemes is not necessary. In the tests reported 
in this paper we have found that in our fully general 
relativisitc implementation, the computational time is 
saved by about 20 %; (3) the density of the unphysi- 
cal atmosphere one needs to build around isolated stars 
when adopting the conservative form of the hydrodynam- 
ics equations can be several orders of magnitude smaller 
than that in HRSC schemes. Associated with this ad- 
vantage the code can be applied for neutron stars with a 
large adiabatic index F = 2.5. 

These facts illustrate that HRC schemes can be useful 
and robust tools for hydrodynamical simulations in full 
GR involving dynamical spacetimes. In addition, their 
suitability over HRSC schemes becomes further appar- 
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